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1. INTRODUCTION 

In previous publications |3IH1) methods to compute the modified Bessel function 
Kia{x) for positive x were described. We complete here this analysis by describ- 
ing analogous methods for the computation of the function Lia{x). With this, 
methods for the reliable computation of a pair of linearly independent numerically 
satisfactory solutions become available which find their implementation in the ac- 
companying paper [Sj. Methods to compute their derivatives are also provided. 

The functions Kia{x) and Lia{x) are solutions of the modified Bessel equation 
for imaginary orders 

x^w" + xw' -h (a^ - x^)w = 0. (1) 



^Present address: Departamento de Matematicas, Estadfstica y Computacion. U. de Cantabria, 
39005-Santander, Spain 
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The function Kia{x) finds appHcation in a number of problems of physics and ap- 
pHed mathematics [7]. The function Lia{x) is a real valued numerically satisfactory 
companion to Kia(x) in the sense described in pp. 154-155. 

In terms of the modified Bessel function of the first kind Ii, {x) , the solutions are 
defined as: 

TT 1 

Kia{x) = r [I-ia{x) ~ ha{x)\ , Lm(x) = - \I -ia{x) + ha{x)\ , (2) 

li smh(7raj I 

with Wronskian W \Kia(x), Lia.{x)\ — \/x. 

Both Kia{x) and Lia{x) are real solutions for real x > and a e M. Because 
they are even functions of a, in the sequel we will consider a > 0, although this 
restriction is not present in the code. 

In Section 2, we describe the different methods of computation considered, namely: 
series expansions, asymptotic expansions for large x, Airy-type uniform asymptotic 
expansions, non-oscillating integral representations (including a discussion on the 
quadrature rule) and a continued fraction method. We avoid duplicating informa- 
tion already given in previous papers; in particular, the references |Z| |Hj provide 
information required for building the algorithms of the accompanying paper j^l . A 
few misprints in are corrected. 

In Section 3, we include a discussion on the dominant asymptotic behaviour of 
the functions. These exponential dominant factors can be taken out, leading to 
scaled functions which can be computed in a much wider range. The algorithm 
described in the accompanying paper ^ offers the possibility of computing scaled 
and unsealed functions. 

2. METHODS OF COMPUTATION 

In [3 IH] methods are described to compute the Kia{x) for different regions in 
the (a;, a) plane. In particular, we considered series expansions ^B], asymptotic 
expansions for large x Eq. 9.7.2), uniform asymptotic expansions for a ~ a: 
(O E] and pg. 425). Also, non-oscillating integral representations ^| Ej 
are available. Similar techniques are available for the computation of Lia{x) and 
the derivatives K[^{x), L[^{x). In addition, a continued fraction method can be 
applied for the computation of Kia (x) and K-^ (x) . Those techniques generally give 
at least two alternative methods for computing the functions in the (x, a) plane for 
moderate values of x and a; therefore, we can always compare different methods 
of computation for testing their accuracy. The selection of one or another method 
of computation in a given region will depend on the range of applicability of each 
method and its efficiency. 

We now describe the different methods of computation which are used in the 
programs. 

2T Series expansions 

Series expansions can be built which properly describe the solutions near the sin- 
gular point (x = 0) of the defining differential equation The idea, as described 
in ^1 [7], is to substitute the Maclaurin series for Ii,{x) ( 1 , Eq. 9.6.10) in Eqs. 
The following expansions are obtained 
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where 



oo oo 
fc=0 fc=0 



Ua{x) = n{a) ^ rfeCfc , L[^{x) = n(a)| ^ 

fe=0 fc=0 



■rk_ 
2 

2/*; 



Cfe 



(3) 



and m 



n{a) 



27Ta 



fk 



sin(0a^fc - Qln(a:/2)) 
(a2(l + a2)...(fc2 + a2))V2^ 



/fc/j^fe = ^ tan((/;>a,fc - a\n{x/2)), with (/>a,fc = arg (r(l + fc + ia)). 



(4) 



(5) 



The coefficients fk and differ from those in 7 by a constant factor (for fixed 
a). The new normahzation shows explicitly (Eqs. l|2Jl) the dominant exponential 
behaviour n(a) and l/n(a) as a ^ 00 (~ e^'^"/^). 

An efficient method to compute the coefficients was described in |71 ^j; this 
method is based on the fact that both fk and satisfy the three-term recurrence 
relation 



+ a^)rk - {2k ~ l)rfe_i + rk-2 = 0. 



(6) 



Perron's theorem is inconclusive with respect to the existence of minimal solutions 
for this recurrence relation; anyhow, the second equation in jsj confirms that neither 
fk nor Tfc are minimal solutions. Therefore, forward recursion will be numerically 
stable. Starting values can be computed taking into account that arg r(l + ia) = 
(To (a), where (To(a) is the Coulomb phase shift, for which Chebyshev expansions are 
available for double precision 3 . Namely, we have: 



rg = cos[CTo(a) — aln(x/2)], 
1 

ri = ^ _^ i {cos[cro(a) — aln(a:/2)] — asin[(To(a) — aln(x/2)]} , 

and 

/o = isin[CTo(a) - aln(a;/2)], 
1 

A = —rr, — 3T {sm[cro(a) - aln(x/2)] + acos[cro(a) - aln(a;/2)]} . 
a(l + a ) 

These formulas correct two misprints in (Eqs. 12 and 13). 
Series can be used for x/a small. See [H], Section 2 . 
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2.2 Asymptotic expansions for large x 

Asymptotic expansions for large x are available from the known asymptotic expan- 
sion ofl^ix) (HI, Eq. 9.7.1): 



ir.(.) = (^)^'"e-^^EM+7. 




(9) 



V2^ \(^; {2xf 

where (ia, m) is the Hankel symbol, which satisfies 

1\2 
k + -\ 



{ia, k+l) = — ^ (ia, k) , (m, 0) = 1. (10) 

fx, ~\~ 1 

Bounds for the error terms (7„, (5„) can be fomid in 11 , Pg. 269, Ex. 13.2. 

As discussed in '7' the numerical performance of the asymptotic expansion for 
Kia{x) is of more restricted applicability than for the case of the evaluation of i^,y(x) 
for real v. Furthermore, the continued fraction method described in 7, covers the 
region where this expansion is of numerical interest. For this reason, the continued 
fraction method is the preferred algorithm for the evaluation of Kia{x) and K[^{x) 
for moderate values of a. On the other hand the asymptotic expansion for Lia{x) 
turns out to be accurate in a wider region, which is a fortunate situation given that 
the continued fraction method is not available in this case. See 0, Section 2, for 
further details. 

Asymptotic expansions for the derivatives are also available by differentiating 
Eqs. 0- 

2.3 Airy-type uniform asymptotic expansions 

The Airy-type asymptotic expansions for Kia{x) can be found in and ^I] pg. 
425; the analogous expansions for Lia{x) and Kl^{x) are also available P], while 
the expansion for L[^{x) can be derived in the same way. We summarize here the 
main features needed for the computation through these expansions, neglecting the 
error terms. Further details can be found in 0^^] and P I18|. 

The expansion for Kia{x) and Lia{x) in terms of Airy functions (Ai(z), Bi(z) 
and their derivatives) reads 



K,,{az) = "V!*^^^^ [Ai(-a2/3c)F,(C) + -^Ai' {-a^^OGaiO 



(11) 



2^ 

where 



Faio - f:i-y^ , G.(c) - f:(-)^Hi^, (12) 



s=0 s=0 
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as a — > oo uniformly with respect to z € [0, oo). Error bounds for the asymptotic 
expansions of the Kia{x) and Lia{x) are given in 0]. 
The quantity C is given by 



2^3/2 ^ log 1 + ^} - - yi^, < z < 1, 

(13) 

2(_()3/2 = 1/^2 - 1 - arccos-, z > 1, 



and 



Of course, it is crucial to compute accurately Eqs. for small (. For this, 

series expansions around z = 1 can be used. 

The evaluation of the coefficients near the turning point z = 1 (which is our region 
of interest) can be performed via Maclaurin series expansions of the quantities i/i, 
Qs and 6s(CHl) in terms of the variable 77 = 2~^^^( (see ^ and jlSj for further 
details) . 

Asymptotic expansions for the derivatives can be found by differentiating Eqs. 
((TT)l . In this way: 



Ki{az)^2^'' 



t/2 



za^^) 



'k^{-a^/H)Pa{C,) + -^Ai(-a2/3C)ga(C)' 
Bi'(-a2/3c:)P,(C) + ^Bi(-a2/3^)Q^(^) 



(15) 



where Pa{0 ^^id Qa{0 be written in terms of -Fa(C)i Ga{C,) and their derivatives 
and they have asymptotic expansions 



PaiO - f:(-)^^ , Qa{C) - f:i~r^, (16) 

^ Qi ^ Qi 

s— s— 

whose coefficients can be obtained from the computed coefficients and hg (in 
Taylor series around C — 0) through the relations: 



c.(C)=a«(C) + x(C)&s-i(C) + &U(C), .^7^ 

rf«(C) = -x(C)«.(C)-<(C)-C^.(C), ^ ^ 

where 

x(C) = 0'(C)/0(C), (18) 

The prime in Eqs. p7|) and denotes the derivative with respect to C- Using 

Eqs. p7|l the coefficients Cs and dg can be computed from the coefficients and 

&s . Details on the evaluation of ag and hg are given in , where an explicit Maple 
algorithm is given for the computation of Og and hg for s = 0, 1, 2, 3. 
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By computing the Wronskian relation for the modified Bessel functions and using 
the Wronskian for Bessel functions, it is easy to derive the relation 

FaiOPaiC) - ^GaiOQaiC) ^ i (19) 

a 

which is a useful relation for checking the correctness of the approximations for the 
coefficients in the asymptotic expansions. 

An algorithm to compute Airy functions of a real variable is needed for the 
computation of these asymptotic expansions. In the routines (H] we use Algorithm 
819 0. 

These Airy-type asymptotic expansion are applied in [2| in a broad region around 
the turning point line a — x. 

2.4 Non-oscillating integral representations 

Paths of steepest descent for integral representations of the modified Bessel func- 
tions of imaginary orders and their derivatives are given in 17 . Apart from their 
application in asymptotics [S], these integrals are useful for building numerically 
stable (non-oscillating) integral representations for Kia (x) and Kl^ {x) , as described 
in We complete here the analysis in by providing analogous expressions for 
the computation of Lia{x) and L[^{x). Additionally, we study further transforma- 
tions of the integrals which enable us to obtain integral expressions suitable for 
computation by means of the trapezoidal rule. 

2.4.1 Monotonic case (x > a). We have the following integral representations in 
the monotonic region Q 



where 



cos 6 



cosh r - 1 + 2sin^ ^(6* - a) 

cos (T 



(20) 



A = a; COS + , a = xsin9, sina = (sin0 ] (21) 

V sinh T / 

and E [0, 7r/2), a E (0, 9]. The dominant exponential term (e^'^) has been factored 
out. The argument of the exponential in the integrand is 

$(r) = (coshr - l)coscr + 2sin (-^—^ ^in (—^^ + (o" - 6*) sin 6*. (22) 

This formula corrects a misprint in 7 (Eq. 33). The difference 9 — a can be 
computed in a stable way for small values of t by using the expression. 



sm 



sm f 



1 - 



cos f 



cos (7 



sinh 



(23) 



sinhr 

together with the definition of a (|21|) and specific algorithms to compute cosh(T) 
and 1 — sinh(r)^/T^ for small r. 
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The non-oscillating integral representations for Lia{x) and its derivative can be 
written after factoring the dominant exponential contribution as: 



e 

2^ 



+ 00 



dr 



(24) 



where 



7(ct) =2sin^^sin^-^ + (cr-6l)sin6' 

■q = 2x[cos e + {e- 7r/2) sin 9] = 2x (^1 - {a/x^ ~ § arccos(f )) 



(25) 



and using Eq. 



da 

—— = tan (T 
dr 



coth T 

T 



(26) 



The first integral is dominant over the second one for large values of the param- 
eters and a/x not too close to a = x. As a — s- a; both integrals become of the same 
order. 

Similarly, we have the following representation for L'^^{x): 



-e+n 



2n 



cos ere' 



''^'^"Ua + (1 - e-2'^")e- 



e-^*M/j(r)dT 



(27) 



where 



w N . fcoshr 1 1 

hir) = sino- — 28) 

T smh T 



These integral representations for Lia {x) and {x) can be used for checking the 
computation of these functions in the monotonic region. They are not used by our 
algorithms |B] because the Airy-type asymptotic expansion f Section I2.3|l and the 
expansion for large x f Section 12. 2|l are sufficiently accurate for this functions and 
they are faster to compute (see (Sj, Section 2). 

2.4.2 Oscillatory case (x < a). The non-oscillating integral representations for 
the oscillatory region are more difficult to evaluate numerically than those for the 
monotonic case. Indeed, as it was discussed in the steepest descent method 
leads to three integrals, which have to be computed separately. However, as we 
later discuss, for moderately large a it will be enough to compute a single integral. 

In 12] , the following formula was obtained: 
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K,a{x) =e-W2 



I COS X + sin X3— I dr 



dr 



+ 



1 r 

ih vra / ^ , 

J /J — tann ; 



da 



■r / I COS Y smh p + sm y cosh p — — dr 

sinhvra / , , \ ^ ^ dr ' 



(29) 



1 

sinh na 



3tt/2 



COS X sinh p - — h sin y cosh p ) da 
da 



where x — ^ sinh /i — a//, cosh fi ~ j:, /i > 0, 



and 



^'(r) = xcosh T cos a + a ( cr — -tt ) ,p(r) = —'!'(''') + an 



(r — /x) cosh p + sinh p 
sinh T 



(30) 



(31) 



Notice that each of the three integrals in Eq. H29() can in principle be integrated 
with respect to any of the two variables a and r, taking into account Eq. H31|l 
together with the fact that the integration path r(cr) is such that r(0) = +00, 
t('k/2) — fi, t{'k) — fj, — tanhp > r(37r/2); however, as discussed in J? there are 
strong numerical reasons for the selections made. In particular, the third integral 
is performed with respect to a (which requires numerical inversion of (|31l) ) to avoid 
the singularity of da/dr at T(37r/2). As explained in ^ the numerical inversion of 
()31|l in the interval a € [tt, 37r/2] can be efficiently performed in parallel with the 
numerical integration. 

Similar integral representations exist for K'^^{x), Lia{x) and L'^^{x). We have: 



_ „-7ra/2 



where 



{cosx A{t) + sin xC(t)) dr 



sinhTra / (cos x cosh p A(t) + sin x sinh pC(t)) dr 

J /J-— tanh fi 
/•37r/2 

sinhTTo / (cosxcoshpB(T(CT)) +sinxsinhpD(T(CT)))dcr 



A{t) = - cosh T cos (7 + sinh t sin , B{t) = A{t) ^ , 

C(r) = - sinh t sin cr - cosh r cos cr^ , D(t) = C(t) ^ . 
In addition, integral representations for Lia (x) and its derivative are: 



(32) 



(33) 
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Lia{x) — ^-jjT 



/2 



l-e 



-*(r) 



sin Y — cos y—r- | dr 



and 



tanh ^ 
37r/2 



sin X sinh p — cos x cosh p — ) dr 

dr 



dr 



sin X sinh p cos y cosh p da 

da 



(34) 



„7ra/2 



1 - e 



-27ra 



e-*(^) (sin X ^(t) - cos x C(t)) dr 



/ (sin X cosh p A(t) — cos X sinh pC(t)) dr 

J tanh 
<.37r/2 

/ (sin X cosh p i3(r(cr)) — COS X sinh p£>(T(cr))) do- 

^ 7r 



(35) 



Notice that the dominant exponential behaviour has been factored for both the 
functions Kia{x) and Lia{x) and their derivatives, which coincides with the ex- 
ponential behaviour of the uniform asymptotic expansion. This is an interest- 
ing feature when computing scaled functions in order to avoid overflows and/or 
underflows in the computation. After factoring the dominant exponential terms 
^•g±7ra/2-j^ the overflow and/or underflow problems are eliminated; notice, how- 
ever, that when computing the integrals over finite intervals we should evalu- 
ate sinh(jo)/e°'^, cosh(p)/e"'^ for Lia{x) and its derivative and sinh(p)/ sinhoTr, 
cosh(p)/ sinhaTT for Kia(x) (and K[^{x)) instead of computing the hyperbolic and 
the exponential separately (otherwise, overflows will take place for moderately large 
a). For this reason it is convenient to use the expressions 



coshp ^ 1 + e sinh p 

sinhTra i — g-^irai sinhvra 



1 - e-^--^ 



(36) 



in the evaluation of Eqs. 129|l and 132|1 and to proceed in the same way for 
e~'^° coshp and e^'^'^sinhp in Eqs. (|34() and 135() . Notice that in the oscillatory 
region p > and that for large a and x both e'^^ and e~^'^" will underflow. These 
underflow problems can be easily avoided by neglecting these exponential terms for 
large parameters. 

In addition, when both exponentials become negligible, the integral over sigma 
becomes negligible and the remaining two integrals can be approximated by only 
one integral. We can write 



10 • Amparo Gil, Javier Segura and Nico IVI. Temme 



Lia{x) 



g-ira/2 


/>oo 






J To 


g-7ra/2 


/>oo 






J TO 


g7ra/2 

27r 


{j 


/•OO 

e 

TO 


g7i-a/2 

27r 


{j 


/•OO 

e 



-*(r) 



da . 

COSY + smx— dr + Cfe 
ar ' 



-7ra/2> 



e-*(^) (cosx^(r) +sinxC(r))dr + C'(e-""/2) 



' sinx-cosx— 1 dT + 0(e-™/2) 
ar / 



(37) 



) {smxAir) - cosxC{t)) dr + Oie-"""/^) 



where 



To = — tanh/i. 



These approximations can be used for moderately large a, which is the region 
where integrals for the oscillatory case are employed in the code 0. 

It is however useful to have the complete expressions for testing the rest of the 
methods. The computation through quadrature using Eqs. (|20f) . H24|) . (|27|l . (|29|l . 
(EH, (EH and provides a way for computing the functions in the whole (x, a) 
plane, except close to a = x, where the integrands become non-smooth. For this 
reason, they have been used for checking the algorithm, although in the oscillatory 
region only Eq. (|37|l is necessary when building the numerical algorithm (Bj . 

2.4.3 Quadrature rule. As reported in Goodwin jlU| . the trapezoidal rule is a 
very efficient method of computation of integrals f{x)dx for rapidly decaying 
integrands f{x); in particular, it is know that the error decays as exp(— (vr/ft,)^) 
for integrals of the type g{x)dx with g analytic in {z e C : \Qz\ < n/h}. 

After appropriate changes of variable, similar arguments follow for integrals over 
finite intervals with a smooth integrand ^| . 

The semi-infinite integrals in this Section are appropriate for their computation 
by using the trapezoidal rule, because they decay as a double exponential as t — > 
+00. On the other hand, the integrals over finite intervals show abrupt variations 
as a — > x, particularly in the oscillatory case, but under an adequate change of 
variables they can be also computed efficiently by means of the trapezoidal rule. 
For finite integrals, we consider a change of variable in order to map the finite 
interval [a, b] into (— oo, -l-oo) and a successive change to improve the convergence 
of the trapezoidal rule | 14l I15| : namely, we consider the following transformation: 



f{x)dx 



, , , , , ,s {h — a) coshi 

g{t)dt, g{t) = f{x{t)) ^ J 
-oo 2 cosh (smh t) 



(38) 



x{t) 



+ tanh(sinhi) 



And the integral is discretized by means of the trapezoidal rule with equal mesh 
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size: 

/+00 +00 
g{t)dt = h g{nh) + e , (39) 

n— — 00 

where the error e is expected to decay very fast as the raesh size is decreased because 
the integrand is analytic and its decay is doubly exponential. We use a trapezoidal 
rule which halves the mesh size until the prescribed precision is reached; the same 
rule controls that the truncation of the infinite series H39|l gives an error well below 
the accuracy claim. 

Regarding the semi-infinite integrals, we use a change of variable to transform 
the integration interval [a, +00) to {—00, 00). We consider the following change of 
variables to perform this map. 

t(2;) =a + sinh-^(e''). (40) 

The additional change x = smht improves the convergence of the trapezoidal rule. 

It is observed that, typically, no more than 8 iterations of the trapezoidal rule 
are needed, which means that the integrands are evaluated at 2^ + 1 = 257 points 
in the worst cases. This is the typical number of iterations for the evaluation 
of /J^j^ dx by means of a recursive trapezoidal rule when double precision 
accuracy is demanded. This fact confirms that the above mentioned changes of 
variable are adequate for the computation of the integrals for the modified Bessel 
functions. 

2.5 Continued fraction method 

As discussed in [7] both Kia {x) and iiT-^ can be computed for moderate a by means 
of a continued fraction method, similar to the corresponding method for Bessel 
functions of real orders (see and ^^Ij PP- 239-240). We refer to [J for a full 
description of this scheme. 

As numerical experiments show, this method is competitive in speed with asymp- 
totic expansions for large x (Section I2.2|l and the range of application is larger. 
Therefore, the continued fraction method substitutes the use of asymptotic expan- 
sions for large x. 

3. RANGE OF COMPUTATION AND SCALED FUNCTIONS 

As described in previous sections, the integral representations which were developed 
indicate that the dominant behaviour for the functions when the parameters are 
large is of exponential type. This means that the computations can only be carried 
for not too large values of a and x in order to avoid overflows/underflows in the 
computation. For instance, from Eqs. (|29|l and (|32|l it is seen that for large a 
(a > x), we have Kia{x) ~ e~°'^/^ and similarly for the derivative, while for Lia{x) 
and its derivative (Eqs. (|34|l and (|35|l ') the asymptotic behaviour is ~ e°-'"l'^ . This 
means that to avoid overflow/underflows in the computation, we must restrict the 
range of a in the oscillatory region to 

a < 21n(10"Af)/7r (41) 

where N is either the inverse of the underflow number (when computing Kia(x) 
or its derivative) or the overflow number (for Lia{x) and its derivative); 10" is a 
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safety factor (in the program, we take n = 8). For processors using the IEEE 
standard in double precision this will approximately limit a to a < 440. On the 
other hand, for the monotonic region {x > a) the integral representations show that 
the dominant exponential behaviour is Kia{x) ~ e~^, Lia{x) ^ e"*"^ where A(a;, a) = 
x{cos9 + 9s'm9), sin6 = a/x {9 e [0, 7r/2]), and similarly for the derivatives. This 
means that, in order to avoid overflows/underflows, the range of computation must 
be restricted to: 



\{a,x) = \J x^ - o? + a arcsin(a/a;) < ln(10"A^) 



(42) 



Figure 1 shows the computable range for 10" ~ lO'^"" (typical value for IEEE 
standard double precision) 

Given that all our expressions have the dominant exponential contributions fac- 
tored out, exponentially scaled functions can be defined which are computable in 
wider ranges. Namely, we define: 

K^a(x) = { K[^{x) 



e^^-^^KlJx) x>a 
e,<^^liK[J^x) x<a 



and 



'^(^''^)L,a(x) x>a 
-"^^^L.aix) x<a 



e--^(^'")i^^(x) x>a 



(43) 



(44) 



500 



400 



300 



200 



100 



100 200 300 400 500 600 700 
X 



Figure 1. Computable range for the evaluation of Kia{x), Lia(x) and their derivatives. 

Note that, as in the rest of the article, we are considering positive a because 
Kia{x) and Lia{x) are even functions of a. Of course, when applying the scaling 
factors for negative a, we should replace a by \a\ in the exponential scaling factors 
of Eqs. and In this way, the scaled functions are also even functions of 

a. 

The definitions in H43|l and (|44|) eliminate exactly the front exponential factor 
in the oscillatory region (a > x) from the series and the Airy type asymptotic 
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expansion and in all the {x, a) plane for the the integral representations. In other 
cases, there remains an exponential factor with soft variation. For example, when 
using Airy-type expansions in the monotonic region (neglecting non-exponential 
factors), we have 

K^a{x) ~ e^-"-/2 = (45) 

where 

A = x{cos0 + (0 — 7r/2) sin^^) = \/ x"^ — a'^ + a(arcsin(a/x) — 7r/2), 

which is small for a; ~ a (6* ~ 7r/2); loss of accuracy in the computation of A for 
a; ~ a can be reduced by expanding A in powers of — 7r/2. 

Similarly, an exponential factor remains when rescaling the asymptotic expan- 
sions and the same happens when applying the continued fraction method. In this 
case, we have for x > a; 

K^aix) ~ e^-^ = e-^ (46) 

where 

A = a;((cos^ — 1) -|- ^sin6) = a; — \J — — a arcsin(a/a;), 

which goes to zero as a/a; ^ (0 ^ 0). Loss of accuracy in the computation of A 
for small {a/x small) can be avoided by expanding A in powers of 6. 
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